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down on computer time, a search was made up over only the two variables 
v,W, Since u, the heading angle, is of relatively small significance in 
the computation of the Hamiltonian. The # n consisted of taking all 
possible combinations of v,w where v ranged over the set of values 3,6, 
9,12,15 and w ranged over the values 200,400,600,800,1000 and comparing 
the Hamiltonian for each set with the one generated in the numerical rou- 
tine. It is possible for the search with this size grid to fail to de- 
tect a set of controls that yield a smaller value for the Hamiltonian. 
However, the time factor is critical and this grid size was felt to fur- 
nish a satisfactory compromise. A search relying on gross computation 
can easily become completely unrealistic in terms of the computer time 
required. 


The Hamiltonian 


— — 

(329) H=: A.V 
has a constant walue on an extremal, with allowance being made for round 
off errors, whenever t does not occur explicitly in H. 

At this point, we might mention that there are two types of wifi. 
tions of the control variables in the classical literature, called weak 
variations and strong variations. [+] Weak variations are variations in 
which the | A | are "small" for each time steif strong variations are 
variations in which í | dut | dt is "small". That is, in weak variations 
only values of control near those used are compared but if strong varia- 
tions are considered, then the new control function may not be "near" the 
one used. 

All methods of determining the routes are methods of variations which 
deform a given path. A path which furnishes a relative minimum, if we 


allow only well variations, may not furnish such a minimum if strong var- 
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latlons are allowed, as will be seen in section 6. 

5. The Numerical Routine. 
The routine for determining the route is given in this section. 
Heuristic Discussion: 


Let us guess a set of values for the parameters hy» h We will 


2° 
then use this set of values to determine the control variables for each 
time t by the minimum principle to determine a route. The terminal point 
thus generated will, in general, differ from the desired one. By chang- 
ing the values of hi, h, appropriately, this terminal point will be 
forced toward the desired point Xp Y 


Mathematical Derivation: 


First we see from (4.3) that 


(5.1) SA = Adh, + A,dh,. 


Note that V = 1, and this in turn implies fV = 0. These facts will be 
used in the following equations. 
Next, from the first of the Euler equations (4.4) by taking the 


total differential we find that 


- v sin ud? + v cos ud + (-Av cos u -/v sin u + E du 
(5.2) 


+(-A sinu - #cos u+ f )dv +f dw= 0, 
uv uw 


if we assume that we can change A,” , u, v, w, at fixed X, Y, Z, t. But 


οἵ 


since 
H = -À v cos u -/ v sin u + f 
uu uu 
H = - À) sin u -J cos u + f 
uv uv 
H = f 
uw uw 


(5.2) becomes 
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ABSTRACT 

In this paper the details of computing an optimum route for a sub- 
marine are studied. 

Typical functions representing the listening devices were used. It 
was found that in some cases several extremals existed and it was neces- 
sary to set up tests for the Legendre and Weierstrass conditions. The 
problem is further complicated by the fact that the optimum control- 
variables may lie on the boundary of the region of allowed values and 
further routines must be adjoined for this. Further, cosners may occur 
and in particular the control may move discontinuously from a boundary 
point to an interior point or vice versa. The routines were made up to 
effect a compromise between the need for accuracy and reasonable com- 


putational time. 
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Introduction 


The purpose of this paper was to develop a numerical routine for 
solving the submarine routing problem; the problem is for the submarine 
to choose a course that minimizes the probability that it will be de- 
tected. Several difficulties arise in the numerical solution and sub- 
routines were included to take care of these. 

In type, the problem is a problem of Bolza in calculus of varia- 
tions. It is complicated by the fact that there may be several routes 
each of which appear to be the solution. They all begin and end at the 
desired points and they all satisfy the Euler equations. It is only 
after routines are incorporated to check other conditions, the condi- 
tions of Legendre and Weierstrass, that it becomes clear whether a 
Particular solution is the desired solution. When these conditions are 
not satisfied this fact must be determined and a routine made up to 
determine the controls to satisfy it. 

In general, the route is generated as an extremal, a solution to 
the Euler equations, though the basic principle is that the control must 
minimize the Hamiltonian. The problem is complicated by the fact that 
the control which furnishes the minimum may lie on the boundary to the 
region, as when the submarine is at maximum depth. Several subroutines 
must be adjoined to treat the case when the control lies on one of the 
bounding faces or edges. 

The existence of a corner introduces further complications; at a 
corner the control is a discontinuous function of time. It is necessary 
to make up a search routine for other values of the control which may 
decrease the Hamiltonian, and it is necessary to compromise between the 


demands of computing time and accuracy in this routine. 


I am indebted to my advisor Professor PF. D. Faulkner for the en- | 

couragement and guidance given me. This is a continuation of the study 

L1] which he began. The principal purpose of this study was to analyze 

the numerous difficulties envisioned and to make up numerical routines 

to resolve these and others that developed, so that a solution could be 

generated automatically on the computer. I also express my thanks to 

Professor W.E. Bleick, whose guidelines were followed in programming 

the numerical routine associated with this problem, and to Mrs. Sally 


B. Kline for help when programming difficulties were encountered. 





1. Statement of the Problem. 
The basic problem studied here is the following. A submarine 
is to make a voyage to point x 


located at point x ul, in a specified 


0° Yo T 


time T. Throughout this voyage the submarine is subjected to enemy de- 
tection devices whose capabilities are assumed known statistically. If 
the submarine has previously gone undetected, let the probability of 
detection in a time interval At be approximately f(x,y,u,v,w,t) At, and 
hence the probability of being detected along the route satisfies the 
equation 

er. 19 dp = (1 - p) f(x,y,u,v,w,t) dt. 
The function f is the best estimate of the enemy's detection capabili- 
ties based on information we have about his listening devices, tests we 
have run on our submarines using comparable devices, the distance in- 
volved in the trip, and other information available to us. 


Equation (1.1) can be simplified by letting 


(1.2) z = —ü— 1n(1 - p) 
which gives 
e. 3) z = f. 


Since we are primarily interested in long routes, and the time to 
change depth, speed, and heading angle is assumed small compared to 
total time, it will be ignored. 

2. Equations. 

Because routes of approximately 2500 miles are of primary interest, 
the flat earth assumption will be used. The equations governing the 
system may then be written as 


x = v cos u 


(2.1) v sin u 


= d 
l 


FG y us SW t). 


Ne 
" 





H = c 


limited range. The constant 6) reflects the possibility of detection by 
some passing ship, for example. 

The thermocline is defined as the depth at which the temperature 
gradient (rate of decrease of temperature with increasing depth) is a 
maximum. In equations (2.3), Wo is the depth of the thermocline, S1 is 
the distance from the terminal point c 2 Y.) of the route to the main 
concentration of enemy listening devices which is represented by the 
point (X55y5); and @ is the angle that the perpendicular to the enemy 
shore line makes with the x-axis. The constants 84 »b45,04,d,, and wg are 
chosen so as to simulate, by the function f, conditions as they exist 
in any given situation. 


The two functions f, and f, are intended to be typical of the 


2 
functions which do describe the enemy's defenses. The function E 
representing the passive enemy defense, tends to decrease as the dis- 
tance from the enemy shoreline increases. This mathematical model of 
the enemy's passive defenses also tends to be insensitive in the region 
of the thermocline. The function £5; which represents the enemy active 
defense, was constructed so as to emphasize the enemy's surface search. 
For this reason, this function decreases as the depth increases. 

If submarine routing as described in this paper were to be made a 
Part of naval operations, generating the function f would be a problem 
for intelligence and engineers. Such things as the sea state and the 
nature of the ocean floor in the region where the submarine is operating 
would then have to be included in the function f. These functions would 
also vary with time, reflecting changing sea state, etc. 

The problem now becomes that of determining the control variables 
as functions of time to effect the desired optimization. That is, we 


want to choose the heading, the speed, and the depth so as to go from 


the initial point (xo Yo) to the terminal point CX Y) with the minimum 
probability of detection. 


By (1.1) the probability of being detected along the route is 


(2.5) p(T) = 1 - exp |- «cn | 
and hence this is the quantity we wish to minimize. 
3. Adjoint Equations. 

Let us consider any route and a neighboring route. The neighbor- 
ing route will be generated by replacing u,v,w on the original route 
by u + fu, v + df, w + dq respectively. This generates first-order 


changes in x,y,z satisfying the differentlal equatlons 


dx = - v sin u du + cos u dv 
(3.1) dy = v cos u du + sin u dv 


d'2 = f du +f dv +f du + f dx + f dy. 
u V W x y 


The notation rx e df 





» etc. as used above will be employed throughout 


the remainder of the paper. 

Let us now introduce three Lagrange multipliers \, ¥,~V which are 
some unspecified function of time t. Multiply equations (3.1) by A,Y, 
Vin that order, add the three equations, and integrate the result over 


the interval (0,T). These operations yield 


T : 
(3.2) Án ENE + v sin u fu - cos u dv +/ (fy - 
v cos u u - sin u fç) + V(éz - £ dx - £ dy 


- f fu - £v - £d | dt. 
u V w 


Separating the terms containing the variations of the state variables 


from those containing the variations of the control variables, we get 


10 





T 
| [Ak «od vd - εκ - ερ] dt 
0 > D 


(3. 3) (|-Av sinu "v cos u * Vf | du + |Acos u 
u 


0 
+) sin u +-f ES + VE dw) dt. 
v wW 
Integrating by parts in (3.3) just the terms on the left involving 
derivatives of state variables, we get 


T ; 
Dd ++ dy + Vz] - f [G + Vf ) 
0 0 = 


(3.4) 


+ dy” EVE) + fz V Jet 


which combined with the right-hand equations of (3.3) yields 


T 
| adx +2 dy +Y η - Los ^ + 
0 0 


T 
(3:5) όν ο” AE) + dz Y dt | [cA v sin u +)v cos u 
0 


ΥΕ; du + (Acos u+Psinu+Vf£ ) dv « V£ dw Jar. 
V W 


Now, choose àÀ,»,¥ in such a manner that they satisfy the differential 


equations 
a = -Vf 
x 
Y = O. 


Equations (3.6) are called the adjoint equations of the variation- 


al equations (3.1). With this choice of 4,* ,v equations (3.1) reduce 
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bo pessat 
ade + dy + Ve = [av sin u -yv cos u 
0 0 
T 
9:7) ve, | du dt + É cos u +Y sinu + Ye | dv dt + 
0 
T 
V £ du üt. 
0 


Note that this formula gives us a relation between the terminal values 
of x,y,z and an integral made up of the variations of the control varia- 
bles du, dv, dw. Note also the important fact that the values of 

dx, dy, dz interior to the interval (0,T) are not needed. 


For convenience in the sections to follow, the vector notation 


= A 
A= |r 
Mv 
(3.8) 
v cos u 
—À 
V =|v sin u 


£ 
will be used. The scalar product of these two vectors 


—< Q — 


(3.9) Hs A-*VW 
is called the Hamiltonian. 


Also for future use, let us define the following three particular 


solutions to the adjoint equations 


> 
T 
X 





“i 
— λ2 

(3.10) A, = 5 = 
EO 
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À 
= 3 T 
az = -( s dt 
0 
ὡς 
1 


4. Conditions for a Minimum. 

In this section, the necessary condition for a minimum will be 
given. 

If a path is to provide the desired minimum, it must first be admis- 
sible. 

Admissibility. One requirement for a route to be admissible is that 
it begin and end at the desired points, i.e., x(0) = 0, y(O) = O and 


x(T) = x_, y(T) = y, for some solution to the differential equations 


T T 

x = y cos u 

(ων) > = v sin u 
z = f. 


A further condition for admissibility is that the depth and speed satisfy 


the inequalities 


0 Š w“ 
(4.1) Ww ΠΚ 
0 < v“ ç 
max 


where w and v depend upon the specific class of submarine under 
max max 
considerätion. 
A set of control variables which are piecewise continuous and sat- 
isfy (4.1) are called allowable. Note that allowability is a local con- 
straint, in terms of time, on the set of control variables. A path which 


satisfies the above constraints, (2.1), (4.1), is admissible. 
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Our problem is to find the one route, among all admissible routes, 
such that z(T) is a minimum. Minimizing z(T) in turn minimizes p(T) and 
this is our objective. 

For a path to furnish a minimum, it must be admissible and also sat- 
isfy the following conditions: 

1. Euler Equations. 

2. Legendre Condition. 

3. Weierstrass Condition. 

4. Envelope Condition. 

The envelope condition was not investigated and will not be treated in 
this paper. The order of the conditions as given above is used since 
this is the usual order in which they will be checked in a problem. 

A point that should be kept in mind throughout is that on a path 
which furnishes the desired minimum, the Hamiltonian must be minimized 
at each value of t. This is the Weierstrass-Pontryagin maximum princi- 
ple with a change of sign. To satisfy this one criterion, the Euler 
equations, the Legendre condition, and the Weierstrass condition must all 
be satisfied. 

The Euler equations, a condition on the first derivatives, require 
that the Hamiltonian have a stationary value. 

In the Legendre condition the second derivatives are checked for a 
minimum. | 2 | The Euler equations and the Legendre conditions are both 
local conditions on the control variables u,v,w. 

Finally, in the Welerstrass condition we compare the Hamiltonian 
for the values of u,v,w used with the Hamiltonian for all allowable values 
of u,v,w, for all values of t. 


Euler Equations. For a minimum, the control variables must be chosen 
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so they minimize 

(3.9) Haz di ww 
as compared with all allowable controls, for all t, O€t «T, for some 
solution A to the adjoint system of equations. Let us consider our solu- 


— 
tion A to the adjoint equations in the form 
— — ἅν 
(4.2) N =h, i fay, We Ac 
ZN ES — 
where h,, h,, h} are arbitrary constants and Λι, ^»: 2 are defined 


in equations (3.10). 


We are now faced with the problem of having three constants in our 
— 


solution Á to the adjoint equations. But we may choose one relation 


among the constants hy > h,,; h,. We chose h, to be 1, and then 


2 3 


(4.3) ANNO. AE ede TS, ee Nes 


This then leaves us with the problem of finding the other two constants 
so that x(T), y(T) will assume the desired terminal values ΧΩ» Vor 
If the values of the control variables lie inside the domain of 


allowed values, then the Euler equations 


D AM 


—À v Sin u +Zv cos u +Yf = O 
ou u u 


(4.4) 3H -4 





Acos u + ⁄“ sin u K = 0 


—— = H =-Yf =O 
aw w W 


are the first necessary conditions for the desired minimum. The Euler 
equations when combined with the adjoint equations (3.6) are referred 

to as the Euler-Lagrange equations. Solutions to the equations of motion 
which also satisfy the Euler-Lagrange equations are called extremals. 


The Euler equations,however,do not insure that we effect the desired 
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minimum for H; they are only necessary conditions. It is then necessary 


to investigate additional conditions, the next being the Legendre condi- 


tiom 

Legendre Condition. This is an investigation of the second-degree 
terms of the Taylor expansion of the Hamiltonian. If we can expand H in 
a series in du, dv, dw, valid in some neighborhood of (0,0,0), the first 
necessary condition for H to be a minimum is that R, = R, = n. = 0, as 


stated in equations (4.4). The second necessary condition is that the 


quadratic form 


H H H du 
uu uv uw 


(4.5) (du dv dw) H H H dv 
vu vv vw 


H H H dw 
(wu wv ww 


be positive semi-definite at least; we hope it will be positive definite. 


This condition, which is known as the Legendre condition, [ 3] may be ex- 


pressed 

H > O 
uu 

H H 
uu uv 

H H > 0 
vu vv 

(4.6) 
H Hes H 


uu uv uw 





0 
Hou Boy H. τ 


Hua Bav Ew 
for the case in which the control variables u,v,w all lie interior to 
the domain of allowed values. 
If one or more of the control variables are on the boundary, the 
conditions are altered accordingly. Consider the case in which w = C— 


for some part of the route. The conditions then become 
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W 
(4.7) Hu B 
> 0 
H H 
vu vv 
whenever w = we ae If w = Wa and us - 0, then we must check to ensure 


that H > 0. Of course if H > O the minimum is not on the boundary at 
that point. 
If v = ete conditions equivalent to those above will be used. 


When both w = w and v=v » the conditions read 
max max 


H < 0 
W 

(4.8) H < 0 
V 

Hiu? 0, 


and if H = 0 on the boundary then we must check How? 0, and similarly 
for H, and Hv 

If the Euler equations are satisfied and if equations (4.6) are sat- 
isfied at some interior point of u,v,w space, then these values yield a 
local minimum; H is smaller at that point than at any other point u,v,w 
in some neighborhood. However, the point may not furnish the minimum 
value for H; it is necessary in theory to compare the value with the value 
for all allowable values of u,v,w. This condition, that the value chosen 
minimize H, is called the Weierstrass condition. 

A comparison must be made between the Hamiltonian for the u,v,w gen- 
erated, H(u,v,w), and for all other allowable combinations of control 
variables, say, EE: zt. de I£ H(u,v,w) > Sa T I» then the * .- 
control variables u,v,w must be replaced by the new set c. vÍ, w, 
which yield the smaller value for the Hamiltonian. 


Unfortunately there is no easy way to check this condition. To cut 
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down on computer time, a search was made up over only the two variables 
v,W, Since u, the heading angle, is of relatively small significance in 
the computation of the Hamiltonian. The search consisted of taking all 
possible combinations of v,w where v ranged over the set of values 3,6, 
9,12,15 and w ranged over the values 200,400,600,800,1000 and comparing 
the Hamiltonian for each set with the one generated in the numerical rou- 
tine. It is possible for the search with this size grid to fail to de- 
tect a set of controls that yield a smaller value for the Hamiltonian. 
However, the time factor is critical and this grid size was felt to fur- 
nish a satisfactory compromise. A search relying on gross computation 
can easily become completely unrealistic in terms of the computer time 
required. 


The Hamiltonian 


(3.9) H= N-V 
has a constant value on an extremal, with allowance being made for round 
off errors, whenever t does not occur explicitly in H. 

At this point, we might mention that there are two types of wits 
tions of the control variables in the classical literature, called weak 
variations and strong variations. [4 ] Weak variations are variations in 
which the | A | are "small" for each time steil strong variations are 
variations in which f. | u: | dt is "small". That is, in weak variations 
only values of control near those used are compared but if strong varia- 
tions are considered, then the new control function may not be "near" the 
one used. 

All methods of determining the routes are methods of variations which 
deform a given path. A path which furnishes a relative minimum, if we 


allow only wel, variations, may not furnish such a minimum if strong var- 


18 


iations are allowed, as will be seen in section 6. 
5. The Numerical Routine. 
The routine for determining the route is given in this section. 
Heuristic Discussion: 
Let us guess a set of values for the parameters h 


h We will 


οσο 
then use this set of values to determine the control variables for each 
time t by the minimum principle to determine a route. The terminal point 
thus generated will, in general, differ from the desired one. By chang- 
ing the values of hi; h, appropriately, this terminal point will be 
forced toward the desired point Xp Y 


Mathematical Derivation: 


First we see from (4.3) that 


(5.1) fA = ^ dh, + A,dh,. 


Note that V = 1, and this in turn implies SV = 0. These facts will be 
used in the following equations. 
Next, from the first of the Euler equations (4.4) by taking the 


total differential we find that 


- v Sin ud^ + v cos ud? + (“Av cos u - v sin u + re) du 
(5.2) 


+(-À sinu - # cos u+* f )ddv +f dw=0, 
uv uw 


if we assume that we can change A,” , u, V, w, at fixed x, Y, Z, t. But 


e, 


since 
H = —À v cos u -/ v sin u + f 
uu uu 
H = =~ À sin u -J cos u + f 
uv uv 
H = Ë 
uw uw 


(5.2) becomes 
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(5:43) -v sin udÀ * v cos ud +H du +H dv + H dwe 
uu uv uw 
Similarly 
cos u dA + sin udy + (-A sinu +/cos u +f ) du 
vu 
(5.4) 
+f (fv +f dy=0 
vv VW 
and since 
H = —ÀAÀ sin u +2Zcos u + f 
vu vu 
H = f 
VV VV 
H z f 
VW TW 
(5.4) becomes 
(5.5) cos udÀ+ sinu f# + H Jdu+H dv+H dy=0. 
vu VV VW 


The third Euler equation gives us 


(5.6) f du +f dv + f dx =O. 
wu wW ww 


But since 


£ = H 

wu wu 

f = H 

WV WV 

f = H 

WW WW 

(5.6) becomes 
5.7) H du«H dv +H dwe 0. 

wu WV WW 


Now consider the equations (5.3), (5.5), (5.7) as three equations 


three unknowns du, dv, dw. 


ν sin udA- v cos u dX 


H du+H W+H dy 
uu uv uw 


(5.8) H du«H dv +H όη Ξξ -cosu dA - sin u d» 
vu vv VW 


H du +H dv +H όν τς 0. 
wu WV WW 
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0. 


in the 





If the determinant of coefficients of du, dv, dw in (5.8) does not vanish, 
we can solve for du, dv, dw by using Cramer's Rule as follows: 


Denote the determinant of coefficients by D, i.e., 


Huu Huy Hw 


(5.9) D= |H H H 


H H H 
| wu WV WW : 
Then 
(v sinud^ - v cos ud") H H 
uv uw 
(- cos u d^) - sin ud») H H 
vv 





(5.10) 


H (v sinudA ~vecosud) H 
uu u 


H (- cos ua -— sinu H 
u V 





D 
Notice that dw could be found in the same manner, but is not needed 


since it is not used in the numerical routine. Next, since from (5.1) 


dA = dh, and d dh, 


du 


-(v sin u dh 


- v cos udh,) (H,, H -H H ) 


1 ww WV VW 
D 





- (cos u dh, + sin u dh,) ο. D - Ην n 


(5.11) 
D 


A 


dv = (v sin u dh, - v cos u dh,) Ou Lar ^ DR p 





D 
* (cos u dh, * Sin u dh,) e lien. s a Hau 
D 
To simplify the form of the equation let us set 
11 | 
ΕΕ | CH H. = Hun Egle + (Ho L1 - 
νν ww WV VW ιν ww 


H v H (- cos υ) / D 


12 
S = | CH H -H H )YC((vcosu +(H H - 
VV WW WV Vw uv WW 
(5.12) HH) (sin v) /D 
s21 = καὶ H H H ) ( ] ) + (H... H 
m. ww Wl aa ee uu "ww 
H ) (cos u) | / D 
wu uw 
22 
S = |a H -H H ) (- v cos u) + (H H - 
vu WW wu VW uu WW 


Hu Hw (sin u) J / D; 


then the above equations may be rewritten 


gu gt! dh, + st? dh, 


21 
dv - 5 dh, i s^? dh, 


(5513) 


We now have the machinery to derive d'x(T), dy CD. In deriving 


—À 
dx(T), we use equations (3.5) and a particular choice for A , namely 


= 1 — 
the Br =10 as defined in equations (3.10). This choice for A in 
0 


equation (3.7) yields 


7 T 
[£ | 2 d'v dt { v sin u Ju dt 
0 


0 0 


22 


which becomes, after noting that d x«o) - 0 and substituting du, dy 


from equations (5.13), 


T 
21 22 11 
(5.14) όχ(τὸ . | [cos u (S dh, + § dh.) - v Sin u (S dh, 
0 
SEE ih )J dt. 
2 
E 0 
Similarly, using A, =| 1}, we find 
0 


1 22 


(5. 159 áy CD [sin u (β΄ dh, + S dh.) + v cos u Cun dh, 


II 
ri 


0 
p 2 
+ S dh.) ] dt. 


Equations (5.14) and (5.15) can be rewritten as 


E 21 š B 
§x(T) | [oos us - v sin u SŠ Jat + 
0 


T 


22 ' 1» 
dh, [cos us - v sin u S jet 
(5.16) 0 


T 
21 11 
dy(T) = dh, [sin us + v Cos u S Jat + 
0 


r 22 12 
dh, [Sin us + v cos u S Jat. 
0 
If we set T 
2] - 11 
All = (cos u S - v sinu S ) dt 
0 


T 
22 à 12 
na | (cos u S - v sin u S ) dt 
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Ass (sin m + v cos u Se dt 


T 
o E 
E 
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12 
A (sin S + voas u S ) dt; 


22. 


i 
El 
ς 


0 
equations (5.16) become 


A dh, 4 Αιρ dh, 


dx(T) 
(5.17) 


Αρ dh, 4 Ayo dh,. 


fy CD 


Equations (5.17) give us the mechanism for a Newton - Raphson iteration 


scheme for correcting hy > ho» if we make the following substitutions 


dx CT Xn - x(T) 


(5.18) 


dy(T) - y(D. 


d 
Note that in equations (5.18) dx(T), dy (T) were equated to the desired 
terminal values minus those which were attained. Using equations (5.18), 


equations (5.17) become 


Xp 7 x(T) 411 dh; + Α12 dh, 


(55 53} 


YT - y CD A dh, + p dh, 


from which we are able to generate corrections for the constants h,, h 


P 


Note above that in equations (5.19) the coefficients A11» Ajo: Αρη» 
A59 are integrals with respect to time. In the numerical routine, this 
integration is accomplished by a Runge-Kutta numerical integration scheme, 
contained within which is a Newton-Raphson iteration to generate the neces- 
sary changes in u,v,w to accomplish the integration. 

The mathematical basis for the Newton-Raphson iteration scheme to 
generate corrections for u, v, w is the following. From the equations 


H du + H dv + H dw = - H 
uu uv uw 


(5.20) H du +H dv+H dw 
vu vv vw 


u 
[ 
am 
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we find 


5:21) 


n 


dv 


uu 


vu 


dw = 


where 





The iteration scheme then has the form 


ú =u 


n n-1 
(5.22) ES UR V 
Wn ^ *h-1 


W 
- H H 
V VW 
D 
H - H 
uv u 
H - Ἡ 
vv V 
H -H 
WV W 
D 
H H 
uv uw 
H H 
vv VW 
Hoy Hw : 
4 du, 
+ dva > 
+ dw 


where n is the index for the Newton-Raphson iteration; du, dv, dw are 


those found in equations (5.21). 
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6. Control Variables on the Boundary. 


In section H, it was noted that for a path to be admissible the con- 


trol variables v, w had to satisfy the inequalities 


|^ 


0 4 


w w 


(4.1) max 


< < 
esas n UU ° 
0 E max 


Let us consider the situation in which the depth w assumes the max- 
imum depth Wise for part or all of the route and the other bounded con- 
trol variable v remains within its prescribed bounds. We must amend the 
routine for determining the controls as follows. We take the maximum 
depth WEE of the type of submarine being considered and read this infor- 
mation into the program of the numerical routine which calculates the 
route. If the Newton-Raphson iterations as established in (5.22) yield 
a value for w which exceeds Wax? Ye set w equal to es. and generate 
u and v, using the following iteration scheme. From our Newton-Raphson 


iteration equations, 


i 
I 
= 


H du + H dv 
(6.1) u uv 

H du + H dv 

vu VV 


ii 
i 
i 


we get corrections to the controls, 


- Ἡ Ἡ 
u uv 
- Ἡ Ἡ 
V VV 
du = 
H H 
uu uv 
un H 
VV 
(ius 2) 
H - Ἡ 
uu u 
Hou -H 
dv z y 
H H 
uu uv 
H H 
vu VV 
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Then the iterations discussed in the previous section take the form 
u =u + du 


(6.3) n  n-l n 


= + " 
Vi, emer dv, 


Thus we may modify our Newton-Raphson iteration scheme for the con- 
trols when the minimum value of H occurs for w on the boundary. The 
modified Newton-Raphson iterations generate successive corrections to u 
and v so as to produce an admissible path. The numerical routine is pro- 
grammed so that it is possible for the ac to assume its maximum for 
some part of the route without remaining fixed at maximum depth after 
once assuming it. In section 7 we will see an optimum path which has 
the control on the boundary for part of the path; the submarine travels 
at maximum depth for a while, then comes up. 

i subroutine for calculating u, v when w - m contained in 
lona: I, part c. The above features can be = by examining either 
the flow chart for subroutine BOUNDW or the subroutine itself which is 
given in part G of Appendix I. i 

Modifications similar to those made for w on the boundary would be 
made if v = v " for some part of the route. One additional change re- 


ma 


quired in this case, which was not necessary when w = Wax? ÍS that dv 


be set equal to zero in the numerical routine whenever v - a and 


LN VE v= Vnax and H, > O imply that the velocity is decreasing or 
moving away from the bound and hence the numerical routine described in 
equations (5.22) for generating corrections to the control variables 
would be used whenever this is the case. It should be noted that setting 
dw equal to zero in the case of = Wax Was not required since dw does 


not enter into the equations for the numerical routine in section 5. 


When v = es, and H; < 0, we get from equations 
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I 
I 
= 


H „du + H, @ dw 


(6.4) 


II 
I 
T 


H du r H dw 


the corrections to the control variables u, w as 
i3 


4. 


- H H 
u u 
- H H 
WW 
du = 5 
H H 
uu uw 
Hoa Hoq 
(6.5) 
H - Ἡ 
uu u 
H -H 
dues | τ- w 
Huu Huw 
Huu S 


Using the du, dw found in equations (6.5), we get equations for correct- 
ing u, w 


(6.6) 
W. = W + dw . 
n n 


ΤΕ both v = v and w = w and it is also true that H < O and 
max max v 


H. < 0, we generate corrections to the control variable u by using the 


equation 
uu u 
from which 
L H, 
(6.8) du  --------- 
H 
uu 


and the iterations to correct u take the form 
(6.9) uM = u + du . 


n n-1 


Both of these situations were encountered in generating Path IV 
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which will be discussed in section 8. The results given there will indi- 
cate that the modifications needed when v - Vaga, SERT κ and w = Wapsi 
do not affect convergence of the Newton-Raphson iteration schemes in the 
numerical routine, because we will find that Path IV is admissible. 

The flow chart for the subroutine BOUNDV, which is used when v = 
Vax? “an be seen in part D of Appendix I. Part B of Appendix I con- 
tains the flow chart of subroutine VUW which takes care of the case where 
e»t). and w z wie The result of setting dv = O when v = ες and 
H, < 0 can be seen by examining the flow chart of the numerical routine 
in part A of Appendix I. The deck listings of BOUNDV,VUW, and the nu- 
merical routine can all be found in part G of Appendix I. 

7. Paths. 

Our computations have established the existence of three extremals. 
Ihe three paths all satisfy the Euler-Lagrange equations as outlined in 
section H. 

These three paths can best be compared by listing the contrasting 
points of the three. The areas in which the greatest difference appeared 
among the three routes are the probability of detection p(T), the depth 
w, and the constants (hy >h,) which yield admissibility. 


The following is a list of the results for each path in the three 


areas just mentioned. 


Path I 
-3 
p(T) . 15711 x 10 
w thermocline = 200 feet 
(hy 5h.) (-.00019, .00024) 
Path II 
-3 
p(T) 72 2725010 


29 


ν 
(hy » ο) 
Path III 
p(T) 
w 
(hy »h,) 


Path I gives us an absolute minimum, 
weak or strong variations. Path II gives 
weak variations are allowed. Path III is 
nish a relative minimum under either weak 


called Path III a worsimax path. 


wW up to 683 feet 


max 


(-.00038, .00048) 


.23060 x 107^? 
525 to 725 feet 


(-.00037, .00046) 


i.e., a minimum under either 
us a relative minimum if only 
an extremal, but does not fur- 


or strong variations. We 


If any additional information concerning any one of the above paths 


is desired, copies of the three paths can be found in Appendix II. 


Given 


there is a printout:. of each path with the coordinates (x,y) denoting the 


submarine's location, the control variables u,v,w, and the probability 


of detection p(t) for each time step. 


Note that in Path II w = w 
max 


for nearly all of the route. 


In cal- 


culating this path, the iteration scheme described in equations (6.1), 


(6.2), (6.3) was used. 
terminal point Gp Yap 
Analysis of Paths. 


4 point out the following facts. 


tions: 
1. admissibility conditions, 
2. Euler equations, 
3. Legendre conditions, and 
4. Weierstrass condition. 
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The fact that Path II converges to the desired 
) substantiates the results given in section 6. 
Checks on the conditions as outlined in section 


Path I satisfied the following condi- 





Path II satisfied 
1. admissibility conditions, 
2. Euler equations, and 
3. Legendre conditions, 
but not the Weierstrass condition. 
Path III satisfied 
1. admissibility conditions, and 
2. Euler equations, 
but not the Legendre conditions nor the Weierstrass condition. 

Hence, by the criterion established in section 4, there is but one 
minimum path, that being Path I. Note that the check of the envelope 
condition was not included in this investigation. 

The above results point out an important fact which is often ignored 
or overlooked; the generation of an extremal, by no means guarantees that 
you have the desired minimum. This emphasizes the need for a check on 
all of the conditions for a minimum at each time step of the numerical 
integration scheme for generating the path. If checking after each time 
Step requires excessive computer time, the checks may be performed at 
some appropriate periodic intervals. 

It can be noted at this time, that if we restrict ourselves to weak 
variations as defined in section à, both Path I and Path II are extremals 
which yield relative minima. In contrast to this, analyzing the paths 
and considering strong variations yields the result noted above, namely, 
that Path I is really the only relative minimum among the three paths. 

After Path II, which does not give a relative minimum, was genera- 
ted, a search, as outlined in the discussion of the Weierstrass condition 
in section 4, was used to determine a new set of control variables u,v,w, 


which would sztí5sfy the Weierstrass condition. The resulting paths 
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turned out to converge to Path I. The subroutine search was used for 
this purpose and is contained in Appendix I. 

A similar search over the grid outlined in section H was performed 
when the Legendre conditions were not met in the case of Path III. The 
set of control variables which gave the smallest value for the Hamilton- 
ian in this search were then used to continue the numerical routine. 
Again the sequence of paths produced by the numerical integration con- 
verged to Path I. 

A path was judged admissible if it came within one-fourth mile of 
the desired endpoint; it was felt that further accuracy was not worth the 
computing time. 

To test the convergence of the numerical routine, on a few paths the 
routine was allowed to continue until no further improvement occurred. 
In each of the three paths above, duplication occurred with accuracy of 
at least one-tenth of a nautical mile. Duplication here means the abil- 
ity of the numerical routine to repeat itself after once converging to 
the desired terminal point. 

It has been noted that a spiral pattern of convergence about the 
desired terminal point is present in the computation of each path. It 
is not clear why this occurs but the following is offered as a possible 
explanation. In the computation, we take H, ts O and vary u,v, 
ο τν, hi,h,, but we assume that x,y have the values they assume on the 
path. For the purposes of explanation, let us examine the equation 


H, = 0, from which we get an equation of the form 


End + Hun ths + H, fu + H, dv + H fs + 


H xX + i dy Ξ0 
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when we take its variations. The last two terms in this expression drop 
out in linear problems and can be shown to be negligible if T is small 

in any case. They introduce considerable complication and extra computa- 
tion and hence were discarded. This omission may be the reason for the 
spiral pattern of convergence. It should be pointed out that the terms 
were omitted in only the correction routine. If the sequence of paths 
converges, there is no related error in the path to which they con- 
verge. 

In generating Path II, a subroutine was used within which the head- 
ing was fixed and a search over the depth and the speed was conducted to 
determine which set of values for these two controls gave the smallest 
value for the Hamiltonian" H. These values were then used in the numeri- 
cal routine to insure a start in the proper direction. This method was 
put into use when it was found that poor initial choices for u,v,w caused 
the Newton-Raphson routine to diverge at the beginning of the route. This 
subroutine SEARCH can be seen in Appendix I, part G/ 

To insure a proper start in computing Path II, the subroutine WORSI 
was used. This subroutine is the same as subroutine BOUNDW described in 
section 6 with one exception, that being that w = Wax is replaced by 
w = 500 or some other intermediate value for the depth. This subroutine 
then fixes w at 500 and computes u and v using the iteration scheme des- 
cribed in equations (6.3). The resulting set of values for u and v are 
then combined with w = 500 to make up the initial guesses for the numer- 
ical routine. This subroutine is also given in Appendix I, part G. 

Few problems were encountered in the generation of Path I, but the 
introduction of conditions to cause the submarine to assume its maximum 


depth Wax OT velocity us aggravated the situation and introduced dif- 
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ficulties to make the subroutines listed above necessary. 
It should be realized that this problem contains some real difficul- 


ties, if approached blindly. With the proper background and forethought, 


most of the difficulties can be anticipated and handled, when encountered, 


by methods such as those described above. 
8. Corner. 

This section contains a discussion of a case in which the path gen- 
erated contains a corner. 

A corner appears when control variables u,v,w whlch minimize the 
Hamiltonian are discontinuous functions of t. For convenience, the nota- 
tion 


_ 
U = 


E <S E 


will be used to denote a set u,v,w of control variables. The conditions 
for a corner are, first, that there exists some point on the route, call 
b 3 p —2 sl 
it ti, where two sets of control variables U and U both minimize the 

; : š —1 
Hamiltonian, H, and second that one set of control variables, say, U , 
gives a smaller value for H when t et. whereas the other set of con- 


UNUS yield a smaller H for t 2 ti: These two conditions can be 


Stated as, first 


(8.1) ud) = H(0%) ΠΠ t = t 


and second 


πο «μῶ΄) for t « bn 
(8.2) 


«21: >? 
H(U ) 2 u(U^) for t > ti 


for some neighborhood of ti. 


At a corner the numerical routine for the corrections should be 


amended by adding terms to handle the changes due to the variation of 
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the corner time ti. 

The changes necessary are outlined in the earlier report “Optimum 
Submarine Routing" [ 14]; Section 6. Even without these correction terms, 
the numerical routine generated an admissible path which contained the 
corner, namely Path IV. 

To construct a mathematical model in which a corner could be ex- 


pected the functions described in section 2 were changed in the following 


manner. We replaced 








TS WM Ww 
So 
o t 
Wo 


with w, equal to five hundred feet. The constants in £1; £5» the equa- 
tions representing the passive and active defenses respectively, were 
chosen in such a way that a corner could be anticipated. 

The model used was one in which the passive defense was dominant 
for the beginning of the route, the two would become equal at approxi- 
mately the middle of the route, and the active search would then dominate 
for the latter part of the route. These facts are apparent when one looks 
at Path IV in Appendix II. 

When we analyze Path IV, we see that the submarine proceeds at ap- 
proximately thermocline depth for the first part of the route and then 
changes to w - Ν a; for the remainder of the route. It can also be noted 


that the speed, v, was considerably less than v - V ees until the corner 
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was encountered, at which time v became equal to Vena Whenever the active 
search completely dominates, the submarine goes as deep and as fast as 
possible. 

A check on the Hamiltonian, H, after each time step shows that H is 
constant, within the accuracy of the routine, for the time steps before 
we reach the corner, but is not constant as we proceed beyond the corner. 
It is not clear why H does not remain constant throughout the route, but 
a possible explanation is that the corner was effectively passed before 
it was found, i.e., the numerical routine failed to detect the corner 
when the conditions for a corner were in fact present. The failure to 
find the corner immediately is a result of the grid size used in the sub- 
routine SEARCH to compare the Hamiltonian for controls U generated by our 
numerical routine with the an s a used in the search. As noted be- 
fore, this grid size was decided upon when a finer grid was found to re- 
quire excessive computer time. Considering that it takes a while for the 
routine to find the corner and it takes the routine a certain amount of 
time to settle down after the corner is found, the fact that admissibil- 
ity was accomplished was thought sufficient to justify omission of the 
corner correction terms. 

Notice that in Path IV v - m with p for a part of the route 
before the corner^and then v - VET and w z eee for the portion of the 
route that comes after the corner. Path IV was generated using the iter- 
ation schemes described in equations (6.4), (6.5), (6.6) when v = v 


max 


and w < w and by using (6.7), (6.8), (6.9) when v = v ad w = w . 
max max max 
The admissibility of Path IV confirms the results stated in section 6 for 


v on the boundary and when both v and w are on the bound of their allow- 


able values. 
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9. Observations. 

This section contains observations which may be helpful to a person 
wishing to continue the study of the submarine routing problem. 

In the Newton-Raphson iterations which occur it may be possible to 
improve both convergence and accuracy as follows. For example, in gen- 


erating the control variables let us make up an error function 


2 2 2 
es H, + e, H, + ey". E 
If each successive iteration does not diminish this error function, then 
we Should diminish the preceeding corrections by a factor of say, two, 
or five. The reason is that the Newton-Raphson iteration moves the var- 
iables in the right direction but may overshoot if the linear terms are 
not dominant. The iteration would be terminated when the above error 
function was less than some preassigned value. The incorporation of such 
routines might well improve convergence, save computing time, and improve 
accuracy. The convergence criterion above is derived from the fact that 
Satisfaction of the Euler equations implies HoH H are all equal to 
zero. Similar conditions could be established for the other Newton- 
Raphson iterations for generating corrections to the parameters h,,h,. 
In the iterations to correct these the established criterion would be 
based upon admissibility of the path. By letting (x,y) represent the 


endpoint of the path generated and (Xps Y) be the desired terminal point, 


we could write the condition as 


2 2 
(x - Χο a < m, 
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10. Summary. 

In this paper the submarine routing problem was studied. Functions 
were chosen that seemed to be typical of the functions representing the 
detection devices, both passive and active. Information that could be 
arrived at only through the use of empirical data such as the sea state, 
for example, was not made a part of these functions. 

Using this mathematical model and determining the path from (X9 »Y9? 
to Gens Y p) in a fixed time T, which minimizes the probability of detec- 
tion, p(T), resulted in the generation of three extremals. Examination 
of the paths using established criterions for a relative minimum lead to 
the following results: one path yielded the desired minimum, a second 
Satisfied all conditions except the Weierstrass condition and the third 
path was just an extremal, satisfying neither the Legendre nor the Weier- 
strass conditions. 

Situations were encountered in which the speed, v, or the depth, w, 
or both v and w were on the boundary of allowable controls. It was found 
that if the control on the bound was set equal to the boundary value and 
corrections generated for the remaining control variables, admissibility 
was accomplished just as when all controls were interior to their allow- 
able ranges. 

With a change in the original model for the active defense, condi- 
tions conducive to a corner were established. The numerical routine 
then generated a path which had a corner and was admissible.  Admissi- 
bility here was accomplished without the use of corrections for the 
corner, and for this reason the corrections were not made a part of the 
numerical routine. 


The numerical routine as described in section 5 was programmed in 


38 


Fortran 1960 for a CDC 1604 computer. Both a flowchart and the program 
for this routine can be found in Appendix I, the flow chart in part A 


and the deck listing for the program in part G. 
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Appendix I 


This appendix contains the fiou charts and the 
deck listings for the program of the numerical routine 
and the subroutines that were used to generate the 
paths described in the text. 

A, Flow chart for the numerical routine. 
[5 TARI] 
[Dimension Y on YVARS, YO,DY,C : AK, TAU,X, Y, Z, UP VT WT | 
pL o NE C T 





f Read : XT,T RDA Ε Bi. 01951. S1,C2, 


D2,X2,Y2,Ht ,E2 (FAC, PAUL, WO,LMÀX, 
HZ, UMAX, VMAX — 






T ο --- 
Q = THETA / 57.29 

COS =»COSF (O) 
SIN = SINF(Q) | 
IX2 = S1 COS + XT | 
| i | 


| E s 5 a 
€ m B s 
‘Print XT,T,THEPA,A1,B1,C1,D1,S1,C2, 
| D2,X2,Y2,H1,H2,FAC,FMUL,WO,LMAX, 
` Wi, WMAX, VHA X '''''. ἵν, —— 
NET RC ο ο δε το. 
| c(1) 2 0.0 Ni =T / XSTEP | 
6G -- Xü-N | 
cG)=00 — STPsT/XM| κα: 
er) = 200 N2 = N1 + 1 
¡XSTEP = T / FAC ` 
kc 
Us ο 
|! V = VAV! 
' W = wO | 
e 
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! mE 





| 
---------------.--,.'.Ν''..''"-Ῥ 
| YVARS(I) = 0.0,I = 1,9] 


NY 


=1 


@— 


Call 
SEARCH 








Qs 


N 
IYC(J) = YVARS(J) + C(I) AK(I - 1,J),J = 1,9] 
| Compute DETER, FUM1 ,FUM2,FUM3| 


ss 


|Cali VUW 
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AE A -ᾱ-- - 





<> 


iDY(3) = COSU m = V SINU 
DY(4) = GOSU S22 = V SINU 
|DY(5) = SINU 521 + V COSU 


~ 


| Sii! 
S12 
S11 





¡DY(6) = SINU $22 + V COSU S12). 


PXC? 
iDY(6 ) COSU 512 


edi Pt των. = FIL + FI2 
DY(8) = V COSU 
DY(9) = = V SINU 


| 


| DY(3) = - V SINU S11 
DY (4) V SINU 512 

V 

ë 


mra e m Se a a 


COSU Sil" 







Me — 
AK(I,J) s STEP DY(J),J s 1,9| 


<> ERR 
E 


— 


IWARS(I) = YVARS(J) + (AK(,d) 
19 AK(3,J) + AK(4,d)) / 6,22 1 
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) + 2 AK(2,9) + 
Y | 


TVAR = TVAR + STEP 
TAU(K) = TVAR 
X(K) = YVARS(8) 
Y(K) = YVARS(9) 







| 
| har = WARS(3) YVARS(6) - YVARS(S) WARS) | 


XNUM1 = (XT = YVARS(8)) YVARS(6) + YVARS(9) YVARS(4) 
IXNUM2 = - YVARS(9) YVARS(3) - (XT - YVARS(8)) YVARS(5) 


fee J. OO eee 
= Hi + FMUL XNUM1 / DET 
H2 z H2 + FHUL XNUM2 / DET 





PROB = 1.0 — EXPF(- 0.0001 Z(N2)) 


NA 
( Print PROB ) 


[STOF] 
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B, Flow chart for subroutine VUW. 


e [START| 
; [I = 1 















| 
Compute PFUMI,PFUMZ « 
FUM3,DETER | 
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Co Flow chart for subroutine BOUNDW. X 


START 


b 
W = WMAX 


NE CS 
[Compute HU,HV,HUU,HVU,HVV| 


M a s — 
DEN = HUU HVV - HVU HVU | 
GUMI = =- HU EW + HWV HWV 
G = + 


UNE 


Beh 
E —— al | 
[DU = GUMI / DEN 
[ον 5 Gui DEN 


| 
| = U + DU 
ο αν 
| 


-w EEE 
Compute HU,HV,HUU,HVU,HVV] 


| 
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D, flow chart for subroutine BOUNDV. 


[SEAT | 


[V = VMAX] 
EE. 
| Compute -— ,HUU , HWU, HWW | 
lcm cc --- ` 
[DEN = HUU HWW - HWU HWU 


GUMI = = HU HWW + HW HWU 
GUM2 = - HW HUU + HU HWU 
Ze | 





Ú = U + DU 
W = VI + DW 


— — 





duci. as BL. 
[Compute HU,HW,HUU,HwU,BWW 


L 





IDEN = EUU EW - HWU HWU 
UMi = = HU EWW + HW HwU 
UM = = HW HUU + HU HWU 


| 


Ea 
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E. Flow chart for subroutine SEARCH. 
START 
Dimension VINT,WINT,H | 


TU cmm NECEM KE ο -αα Oo o A E 
HAM = XLAM V COSU + XMU V SINU + FIl + FI2] 





NN 


— ——— o AS D 


VINT(O) = 0.0 
WINT(O) = 0.0 
oe A 


iW zi 


E 


— — — .  — — ο 


Ζ 
[| VINT(M) = VINT(M - 1) + 3.0] 
| 


x, 


5 un 


 [WINT(N) = WINT(N - 1) * 200.0] 





ποπ κ πω 08 mern: — "NNNM — 
IEQLN) 2 XLÀM VINT(K) COSU + XMU VINT(M) SIND + FIL + FI2| 


— N 

IHAMSRCH = B(1,1) 
VSRCH = VINT(1) 
WSRCH = WINT(1) 


48 li 





HAMSRCH = H(M,N) 
>VSRCH = VINT(M) 
|WSRCH z WINT(N 
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F, Flow chart for subroutine WORSI, 


START 





— — 
DEN = HUU HVV =- HVU HVU 
GUMI z - EU EVV + HV HVU 
GUM2 = - HV HUU + HU HVU 














DEN 
DEN 






DU = GUM1 / 
DV = GUM2 









DEN = HUU HVV - HVU HVU 
GUM = =- HU HWV + HV HVU 
GUM2 = = HV HUU + HU HVU 








(Ο ΩΩ 


E 


{ο Fortran Program: Printout of deck of punch cards 
submitted to computer, 


The main numerical routine and various subroutines 
are given here. The equations in the numerical routine 
and the first four subroutines as listed here are the 
ones which gave a corner. The equations in subroutine 
WORST are the ones used in generating the three extremals. 


PROGRAM SUBROUTE 


YVARS(i)=LAM3 YVARS(2)=MU3 YVARS(3)=A11 YVARS(4)=A12 YVARS(5)zA21 
YVARS16)=A22 YVARS(7)= Z YVARS(8)= X YVARS(9)= Y 
XLAM=H1I+LAM3  XMU=H2+MU3 


DIMENSION YVARS(9)9YC(9)»DY[(9)9C(4)»AK(4»9)9»TA4U(400)9X(900)» 
- ¥(90 )52Z(400) 5VT(400) sUT( 400) »sWT (400) 

READ lLoXToTsTHETAsAl 981 2Cl 2D15S19C22929X22Y2 sH1 sH2sFACsFMUL » 
+WOoLMAX oW 1 o WMAX s VMAX 
1 EORMAT (75,0. 2F 4,0. F2.90»F342»F2el15sF5e 4S F4.O0sFTeG6sF5oe 4s F6oe1l»F 4.41 
+ F6.4>F6.0>F4,05F4.25>F4.0512/F4.0>F5.0>F3.0) 

QzTHE]A/ZOTQOIOIENT T951 

COS= COSA LOS 

SIN=SINF(Q) 

X2=S51*CO5+XT r 

Y2=S1%*SIN 

VAV=XT/T 

PRINT 2 
2 FORMAT (C1HO2HXT4X1HT3X5HTHETA2X2HA12X2HB12X2HC12X2HD12X2HS 12X4HWMAX 
+1X4HVMAX) 

PRINT 39XTosTsTHETAsA19B1l9C15Di 5S19WMAX 5 VMAX 
3 FORMAT (3F5e0>F5e1»F4.25>F4el»F504>F5.0,F5.0>F3.0/) 

PRINT 4 
4 FORMAT( 3H C24X2HD24X 2HX24X2HY25X2HH1 7X2HH23X BHF AC2X 4HLMAX 
+3X4HFMUL2X2HW05X2HW1) 

PRINT 5>C2>9D29X29Y29H1>H2»FAC3LMAX >»FMUL9WO > Wl 
5 FORMAT (F6e55F5.1»F6e1»F4.1»2F94 BS F6609 I3»F8e2»*F6605*F640/) 


DEPUNE REUNGE-KUTTA CONSTANTSe 


C(1)=0.0 


C 


C 


C(2)20,5 
C(3)20,5 
C(4)=100 
X(1)20,0 
Y(1)20,0 
TAU(1)20,0 


COMPUTE TIMO STEP LENGTH. 
XSTEP=T/FAC 
NlLET/ASTER 
XN1=N1 
STEP=T/XN1 
N2=N1+1 

GUESS INITIAL VALUES FOR Us Vs AND Ww BEFORE FIRST ITERATION 
U=-.9 ° 
VzVAV 
W=WO 

ENTER LOOPMEFOR GENERATING CORRESTIONSADO AL. HZ. 
DO 14 L=1»LMAX 
TVAR=0.0 

GUESS INITIAL VALUES FOR UsVsWs AFTER FIRST ITERATION 
IF (L-1) 1985,199,198 

198 U=UZERO 
V=VZERO 
W=WZERO 

199 UT(1)=U 
VT(1)=V 
WT(1)=W 
DO 6 I=1»9 

6 YVARS(I)20,0 

ENTER LOOP FOR COMPUTING THE PATH FOR EACH TIME STEP, 
DO 11 K-2»N2 
IF(K-3) 901,902,902 

O02 CA ee SEAR GHA s 51,C1 51, FI1 FI O ο. το: »V»WsSINUsSWO, 

+GlsG2ZsWlsCOSUMQ) 
ENTER LOOP FOR RUNGE-KUTTA INTEGRATION. 
901 DO 88 I=1>4 


f 
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moO S l O 

ΕΤ Y VAR Ito +C (T)%AKUI-—Ls2) 

AL AMP E e C 1) 

XMU=H2+YC(2 ) 

XD=(XT—YC(8))%COS+S1—-YC(9)*%*SIN 

FI1=EXPF(—(—(XD*LOGF (2.)Z/1000.,)) 

FI2=C2+D2*EXPF (-((YC(8)-X2)**2+(YC(9)-Y2)**2)*LOGF(2.)/10000008). 

GCOsSU=COSF(U} 

SINU=SINF(U) 

COSUMQ=COSF (U-Q) 

SINUMQZSINF (U-Q) 

Gi=C1+D1* (W-WO) *(W-WO) 

G2=1.2.t+D1* (W-WO) *(W-WO) 

G3=(+.01)*FI1*(A1+B1*V%V) 

G&2(-Gl1*(2.*D1*(W-WO))4G2*(2.,*D1*(W—-WO)))/(G2*G2) 

G5z(Cl-4D1*(W-WO)*(W-WO))/(1e-*D1*(W-WO)*(W-WO)) 

G7-2letW*W/((Wl*Wl) 

Pl=((1.+D1*%* (W-W0)* (W-W0))%*20*D]1+4.*D1* (W-WO)*D1I*(W-WO))/(G2*G2) 

P2=(-(8B8Bo*D1%(W-WO)*%*(W-WO))%D1)/(G2%*G2) 

P3= (-(CleDl1*(W-WO)*(W-WO)53*2,*D1-4,.*D1*(W-WO)*D1*(W-WO))/(G2*G2) 
42 (G1*8,*D1*(W-WO)*Dl1X(W-WO))/(G2*G2*G2) 

P9zFI2*(-2./(W1*W1*G7*GT7)--8.*W*W/(WlX*Wl*W1*W1*GT7*GT7T*GT7)) 

HV=XLAMX*COSU+XMU*SINU+ (.01)*F]1]1%*(2.*B1*V)*G5*(l.-oe25*COSUMO) 

HU--XLAM*V*SINUTXMU*VX*COSU-G3*G5*(4.25*SINUMQ) 

HVV-2(.01)*2,*B]*FI1*G65*(1e—e25*COSUMQ) 

HVUZ-XLAM?*SINU-XMU*COSU-G5*,01*FI1*2,*B1*V*(,25*SINUMQ) 

HUUZ-XLAM*V*COSU-XMU*VX*XSINU-G3XG5*(,25*COSUMQ) 

HW=G3%*([1l.-.25* (1 COSUMO ))%*G4+FI2X(-20* (W/W1)*(12./W1))/(67*G7) 

HWUZG3*(4,25*(SINUMQ) ) * G4 

HWV=G4*(.01)x*F11*(2.*B1*V)*(1l.-.25*(COSUMQ)) 

HwW- G$*((15.-.25* SOS Uum sai sa!) 

PD5zHVV*HWW-HWV*HWV 

Pó= (-HV)*HWW+HW*HWV 

P7=HVV#*(-HW)+HWVFHV 

DETERZHUUXPS5-HVU* ( HVUX*XHWW- πα Hw eei Vn 

FUM1z2(-HU)*P5-HVU*P6-c-HWU* ( ( -CHV ) *HWV--HW*HVV) 
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* 
= —— en n tr 


cUM2=HUU*P6+HU% (HVU*HWW—HWUXHWV ) +HWU+* (HVU*(—HW)+HWUZ<HV) 
-UMA-HUUXP 7-HVU* CHVUX ( HY ) -HWUHV ) -HU* HVUXHWV-HWUXHVV ) 
CALL SUBROUTINE FOR CORRECTIONS TO U»VsWe 
CALL VUWLDETER»FUMI1»FUM2sFUM3SXLAM»XMU,FI1»FIZsAl»B1»C1»D1»WOs 
+Vl1 sWMAX»VMAX sUsV>W>Q) 
FIV=COEFF OF FIl F2VsCOEFF OF FI2 
IF(K-2) 19592979195 
297 IF (i-1) 19591969195 
196 UZERO=U 
VZERO=V 
WZERO=W 
195 XD=(XT-YC(8))*¥COS+S1-YC(9)*SIN 
FIL=EXPF(-(XD*LOGF (22)/1000¢)) 
F122C24D2*EXPF (-(LYCU8 )-X2) **2- (YCU 9) -Y2) **2) XLOGF (24 )/1000004,) 7 
G1=C1+D1*(W-W0) *(W-WO ) 
G2=10++D1* (W-W0)*(W-WO ) 
G3=(+.01)*FI1*(Al+B1I*VV) 
G4= (-G1%(2.*D1% (W-WO0) )+G2* (2. *D1*(W-W0)))/(G2*G2) 
G52(CI4D1*(W-WO)*(W-WO) )/ C14 *D1* (47WO)* (W-WO)) 
G721e*(W/W1)**2 
COSU=COSF(U) 
SINU=SINF(U) 
SINUMQ=S INF (U-Q) 
COSUMQ=COSF (U-Q) 
Pi=((1etD1¥ (W-WO)* (W-WO) )*20*D14+4e*D1* (W-WO) *D1¥*(W-WO) )/(G2*G2) 
P2-(-(Be*D1*(W-WO)*(W-WO))*D1)/(62*G2) 
P3= (C(CI4DI* (W-WO 3 C47WO) )2,*D174 ,*DY (W-WO) YD1*(w-WO))/(G2^G 
P4z (61*8, *D1* (W-WOY*D1*(W-WO))/(G62*62*G2) 
P9-FI2*(-26/(W1*W1*GT*GT7)-8e Cw*W/ (W1*W1*WI*Wl1*GT7X*GT*GTI)) 
HVEXLAM*COSU*XMUXSINU* (401) *F 11€ (24 *XB1*V )*G5*(1, 74 25 XCOSUMQ) 
HUs-XLAM*VX*SINU-XMU*VXCOSU*G3*G5* ( 425 * SINUMQ) 
HVV2(.,01)*2, *B1*FI1*G5X(1e-625*COSUMQ) 
HVUZ-XLAM*SINUCXMU*COSU-*G5*,01X*FI1*2.,*B1*V*(4,25*5 INUMQ) 
HUUs-XLAM*V*COSU-XMU3* VXS INU*633 65** ( 4 25* COSUMQ) 
HW=G63%(l.-.25%*(COSUMQ) )*G4+F12*(-20*(W/W1)*(102/W1))/(67*G7) 
HWU=G3%(.25%* (SINUMQ))*G4 








HWV=G4*(.01)*FI11*(2.*81*VY)*(1,-.25* (COSUMO)) 

Ην G3*((1le—-e25*COSUMQ) * (P14+4P2+P3+P4) )+P9 

FIV=(A1+561*VxV)*(.01)*(l.-e25*COSUMQ)*G1/G2 

F2VW=l./(l.e+-(W*W/(Wi*W1))) 

FY=FIV*FI1XSINF$LOGF(2.)/1000»+*F2V*(-D2%*EXPF(=( (YC18)-=X2)%*2+ 
πι Υτ (ο  -Υ2) ΔΡ ΟΟΓ ορ ὐυ 12000000009 £ (2. * CYCUO Y —-Y2) *EOGREU2. )7100000., 
+)) 

FXsFYiV*FIlX*SIN*LOGF(2$.)/1000e4F2V*(-D2*XEXPFC-CUYC(8) -X2) **24 
-£YC(3)-Y2)**2) *LOGF(24. )/1000009, )*(2.*(YCCGG8) 7X2) *LOGF(24,)/100000, 
+)) 

DENOM=HVU* (RVU*HWW-HWU*%HWV ) -HUU* (HVV %HWW- HWV*HWV ) +HWU* ( HVV*HWU-— 
rHWVFHVU) 

S11=((HVVXHWW-HWV*HWV)%(-VXSTINU)+(HVUX*HWW-HWV*HWU ) * ( CCOSU)) /DENOM 

S12= ((HVV+HWW-—HWVXHWV ) *(VXCOSU ) -CHVU*XHWW-HWV *HWU) *( -CSINU) )/ DENOM 

S212CCHVUXHWW—-HWUSHWV) XCV*SINU)-c (HUUXHWW-— HWU*HWU ) *COSU) /DENOM 

$S22z2(0IOHVUXHWW-HWU*HWV)*(CC-V) *COSU) Y LHUUXHWW-HWU*X*HWU) *SINU)/DENOM 

DY(1)=—FX 

DY(2)2-FY 

IF(V-VMAX)118,119,119 

118 DY(3)=C0S5U*S21-V*SINUXSI11 

DY(4)2COSU*S22-V*SINU*S12 

DY (5)=SINU*S21+V*C05U*S11 

DY (65)=SINU*X*S22+VX*COSU*SI2 


OQ TO 120 
119° V=VMAX 
. ΡΥ (3) -V*SINUX*S11 
x. py (4)= -V*SINUXS12 
DY(5)= V*COSU*S11 
BY(6)= V¥COSU*S12 


mao DY( /)=FII*FIV+FI2*F2V 
DY (8)=V*¥COSU 
DY (9)=V%SINU 
DO 8 J=1,s9 
S ΡΟ ΞΟΤΕΡΣΟΥΠΟ) 
88 CONTINUE 
DO 9 J=159 


e tr 4 


— —— = 
Qe 
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ld 


YVARSCJ) YVARSCJ) c CAK( 132) 29 XAK (2$ 2) 29 XAK( 34 2) AK (4,9 ) ) 765 


TVAR=TVAR+STEP 
TAU(K)=TVAR 
X(K)=YVARS(8) 
Y(K)=YVARS(9) 
Z(K)=YVARS(7) 
UT(K) =U 
VT(K)=V 
μτιΚ)-ν 
CONTINUE 


DET=YVARS(3)*YVARS(6)-YVARS(5D)X*YVARSI(4) 
XNUM1=(XT-YVARS(8) )*YVARS(6)+YVARS (9) ¥YVARS (4) 
XNUM2=~-YVARS (9) *YVARS (3 )-(XT-YVARS (8) )¥YVARS (5) 


H1=H1+ FMUL*XNUM1/DET 
H2=H2+ FMUL*XNUM2 /DET 
NA 


12 
13 
19 
‚20 

14 


21 


FORMAT (2X2HN25X2HH17X2HH25X5HX(N2)4X5HY CN2) 6X5HZ (N2)) 
PRINT 135»N25H1»H2»X(N2) 3$YCN2) 4ZCN2) 
FORMAT (4245 279,5,279014E13,57) 


PRINT 19 


FORMAT (1HO3X 3HTAU6X LHX7X LHY 7X LHUSX1LHV5X LHW1OX1HZ/) 
PRINT 20e¢(TAU(T)oX(1)5Y(1I) sUT( I) sVICI) Ντ 1) sZt 4 ) 5 I=1yN2) 


FORMATE BE ZIP BEL +3 FE, 1, E 


CONTINUE 


PROB=1;0=EXPF(=0.0001XZ(N2)) 


PRINTZ IS PROB 


FORMAT (1HO2X5HPROB=E13.5/) 


SIOR 
END 


SUBROUTINE VUW( DETER »sFUM1 » FUM2 »FUM39XLAMsXMU sFI1sFI29A1sBl9ClsDls 


+WOoWloWMAXs VMAX s U», VeWsQ) 
DO 5 I=15s7 
DUZFUMI/DETER 


s sas T e 
: E Ta po es É 


name es cec eee 


A pA penan pap. p. 
faena s. 


—sa 


e 


w sire teen 





00 
ET 
. 3000 


4900 
1000 


DV=FUM2/DETER 

DW=FUM3/DETER 

U=U+DU 

V=V+DV 

W=W+DW 

IF(WMAX —W) 39023909392 

IF(VMAX-VW) 1000910001001 

IF(VMAX-V) 3000300093001 

CALL BOUNDV(UsVseWsXLAMsXMUsFI1,A19619sC1lsD1lsWOsFI29Qs VMAX sHVsW1 ) 
IF(HV-0.0) 4000940003001 

RETURN 

ο ΘΟ ΞΤ. 7 

V=VMAX 

W=WMAX 

Οσο 

SINU=SINF(U) 

COSLIMO= COSF (U=Q 

SINUMQ=SINF (U-Q) 

GisCle-Dl1*(W-WO)*(W-WO) 

G2=10+D1*%*(W-WO0)*(W-WO ) 

G3=(+.01)*FI1*(A1+B1*V*%V) 

G4=(-G1*(2.*D1* (W-W0) )+62%* (20X*D1*(W- -WO)))/(G2*G2) 
G5z(Ci«-Dl1*(W-WO)*(W-WO))/(14-*0D1*(W-WO)? (W—-WO)) 
HUZ-XLAM*V*SINUTXMU*V*COSU-*G3XG5*(.25*S INUMQ) 
HUUZ-XLAM*VX*COSU-XMU*V*SINUTG3*G5*(4.25*COSUMQ) 
DI=-HU/ ĦUU 

U=U+DU 

CONTINUE 

RETURN = A 
CALL BOUNDW(UsVoWsXLAMsXMUsFI1sAlsBlsClsDlsWOsFI29W1»sQsWMAX s HW) 
IF(HW-0.0) 5000500093001 

RETURN 

SOSU-COSF OU) 

SINUZSINF (U) 

COSUMQzCOSF(U-Q) 

SINUMQzZSINF(U-Q) 


no sa αἴ, τ. , 
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G1=C1+D1*(W—WO)*(W-WO). 
G2=1,+D1%(W—WO)*(W—WwO) 

G32(,01)*FI1X*(A1481*VX*V) | 

G4z2(-G1*(26 *D1* (W—WO) ) 3 G2* (20 *D1X* (W—WO)))/(G2*G2) 
G5=(C1+D1*(W-W0)*(W-W0))/(la +OTA WO IO) 
G7=letWEW/EWLEWL) 

Piz=((1Le4+D1¥(W-WO) #(W-WO) )*20%D144e¥D1*( WR WO) *D1#(W- WO))/(G2*G2) 
P2=(-(86*D1*(W-WO) ¥(W-WO))*D1)7(G2*G2) 

P3= (-(C1+D1#(W-WO)*(W-WO) )¥2,.%D1—-4e*D1*(W-WO) *D1*(W-WO))/(G2*G2). 
P4= (G1*8.*Dl1*(W-WOD)*D1X(W-WO))/(G2*G2*G2) 
PosFI2*(-2./(W1*W1*GT7*GT7)*8. *wxW/(CW1*W1X*XW1*W1*G7*G7*G7)) 
HV=XLAM*COSU+XMU*SINU+(.01)%*FI1*(2.*B1*V)*G5*(l.-.25*%*COSUMO) 
HUz-XLAM*VX*SINUC-XMU*VXCOSU«-G3*6G5*(4.253*S INUMQ) 
HVVz(.,01)*2,*B1X*FI1*G5*(1«—6525*COSUMQ) 
HVUZ-XLAM*SINU-XMU*COSU-G5*,0]*FIl*2,*Bl*VX*(.,25*SINUMQ) 
HUUzZ-XLAM*VXCOSU-XMU*V*SINU-G34G5*(,25*COSUMQ) 
HW2zG3*(1es-e25*(COSUMQ) ) *G&44FI24(-224,* (W/W1) C19 7/W1))7(G7*GT) 
HWU=G3%(.25%(SINUMQ))%G4 
HWV=G4X*(.01)*FI1*(2.*B1%*V)*(l.-.25%*(COSUMOQ)) 

HWW=  G3*((1.-.25%*COSUMQ)*(P1+P2+P3+P4))+P9 

P5zHVVXHWW-HWVHWV 

P6z(-HV)*HWW-4HWAHWV 

P72HVV* (HW) +HWV*HV l 
DETER=HUU*P5-HVU% ( HVU*HWW-HWU*HWV ) +HWU* ( HVUX*HWV=HWUXHVV ) | 
FUM]=(-HU)*P5-HVUXP6+HWU*( (-HV ) #HWV+HW*HVV) 

FUM2 =HUU*P6+HU* ( HVUXHWW-HWUXHWV ) +HWUX%(HVUX* (—HW) +HWUXHV ) 
FUM3zHUU*P 7-HVU* ( HVU* ( CHW) -HWU*XHV ) -HU* ( HVUX*HWV-HWUX*HVV ) 
CONTINUE 

RETURN © 
END | 
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SUSROUTINE BOUNDW(UsSVsSW*XLAM» XMUSFI1»A15815» C1,» D1; WOs FI2»W1sQ*sWMAX » 
+HW) 

W=WMAX 

COSU=COSF (U) 

SINU=SINF(U) 

COSUMQ=COSF (U-Q ) 

SINUMQ=S INF (U-Q) 

G1=C1+D1*(W-WO) *(W-WO) 

G221e-Dl1*(W-WO)*(W-WO) 

G3=(,01)*FI1*(A1+B1*V¥V) 

G4á-2(l-G1*(2. *D1* (W-WO) ) 62* (2e *D1* (W—WO)))/(G2X*G2) 
G5=(C1+D1*(W-WO)*(W-WO) )/(1e4+D1*(W-WO) *(W—WO) ) 
G7=let+(W/W1)**2 
HV=XLAMX*COSU+XMU*SINU+(.01)*F]1]1*(2.%*B1*V)*G5*(1.-.25*COSUMO) 
HUZ-XLAM*V*SINUT-XMU*VXCOSU*G3*65(.,25*S INUMQ) 
HVV2(.,01)*X2,X*B1X*FI1XG5*(14.—.25*COSUMQ) 
HVUZ-XLAM*SINU-XMU*XCOSU-G5*,01X*FI1*2,.*Bl1*V*(,.25*SINUMQ) 
HUUZ-XLAM*VXCOSU-XMU*VX*SINU-G3XG53(,253 COSUMQ) 
HWzG3*(1.-.425*(COSUMQ) ) *GA&«-FI2* (C2, * W/W1) *C 1, /W1) ) 7 CGT*G7 ) 
DEN=HVU*HVU—HUU*HVV 

GUM1=-HV*HVU+HU*HVV 

GUM2zHVU* (HU) XHV*HUU 

581 I=1,7 

DU=GUM1/DEN 

DV=GUM2/DEN 

U=U+DU 

V=V+DV 

EOSU-ZCOSF (U) 

SINU=SINF(U) 

COSUMQ=COSF (U-Q) 

SINUMQ=SINF (U-Q) _ 
Gl2Cl«D1*(W-WO)*(W-WO) l 
G2=1e+D1*(W-WO) *(W-WO) 

G3=(.01)*F11*(A1+B1%*VxV) 
G4=(-G1*¥(2e*D1* (W-WO) )+G2*(2e*D1*(W-WO) ))/(G2*G2) 
G5=(C14+D1*(W-WO)*(W-WO) )/(1e+D1*(W-WO) * (W-WO) ) 
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ee 


ο, 
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Vai re REM. s ` dee Mess gas š a 
ο απ οσο κο ασ νο ο a ANDINO Aia doe NG 





π᾿ 


+Wi) 


HV2XLAM*COSU-XMU*XSINUT(,O1)*FI1*(24.*81*V)*G5*(1, —.,25 *COSUMQ) 
HUz-XLAM*VX*XSINU-XMU*V*COSU-G3XG5*(4.25*SINUMQ) 
HVV=(.01)*2.*B1*F11*G65*(1l+.-.25* COSUMO ) l 
HVU==XLAM*S INU+XMUXCOSU+G5O*.01%*FI1*2.*B1*V*(.2595*SINUMO) 
HUUZ-XLAM*V*COSU-XMU*V*SINUTG3*XG5*(4,25*COSUMQ) 
DEN=HVUFHVU-HUUFHVV 

GUM1l-2-HVXHVU-HU*HVV 

GUM2zHVU* (HU) -*HV*HUU 

CONTINUE 

RETURN 

END 





SUBROUTINE BOUNDV(USVsS Ws» XLAM»s XMUSFII1sAl1»5Bl1»Cls»D1s»sWOsFI2»Q» VMAX3 


V=VMAX 
COSU=COSF(U) 

SINU=SINF(U) 

COSUMQ=COSF (U-Q) 

SINUMQ=S INF (U-Q) 

G1=C1+D1% (W-WO)*(W-W0) 

G2=1letD1*(W-WO) *(W-WO) 

G3=(.01)*FI11x(A1+81*V*V) 

G4=(-G1%(2.*D1%* (W-W0))+G62*(2.*D1%*(W-W0)))/(G2*G2) 
65=(C1FDIFUW=WO 1% (WO 9/11 +DIFIW-WOYE UWZWO)) 
GT=1e+WAW/(W1XW1) 

Pl=((1.+D1%* (W-W0)% (W-W0))%2.*D1+4se*D1* (W-WO) *D1*(W-W0))/(G2*G2) 
P22(-(8. *D1*(V-WO)*(W-WO))*D1)/(62*G2) 

P32 (-(Cl14D1*(W-WO)*(W-WO))*2,*D1—44.*D1*(W-WO)*D1*(W-WO))/(G2*Gj 
P4= (Gl*8.*Dl*(W-WO)*D1*(W-WO))/(G2*G2*G2) 

PO=F12%* (-20/(W1*W1%*G7*G7)+8.*W*W/(W1*+W1*W1*W1*G7*G7*G7)) 
HUZ-XLAM*VXSINU-XMUX*VXCOSU-TG33465*(,25*S INUMQ) 

HUUZ-XLAM*VX COSU-XMU*VXSINU-G3:465*(4.25*COSUMQ) 

HW2G3*(1.- 425% (COSUMG) 18 GH E12 CORE EE 
HWUzG3*(4.25*(SINUMQ) ) *G4 
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EWw=  G3*((l.-.25%*COSUMO)*(Pl1+P2+P3+P4))+P9 
DENSZHWU*HWU-HUU*HWW 

GUM]1=-HW*HWU+HUx*XHWW 

GUM2=-HU*HWU+HW *HUU 

DO 981 I-1s7 

DU=GUM1/DEN 

DW=GUM2/DEN 

U=U+DU 

W=W+DW 

COSU=COSF(U) 

SINU=SINF(U) 

COSUMQ=COSF (U-Q) 

SINUMQ=SINF (U-Q) 

G1=C1+D1*(W-W0)*(W-WO) 

G221.-D1*(W-WO)*(W-WO). 

G3=(+.01)*FI1*(A1+B1*%*VX%V) 

G4&z(-G1*(2s. *Di*(W-WO))-74G2*(25s*D1*(W—-WO)))/(6G2*G2) 
G5=(C1+D1%*(W-WO0)*(W-W0))/(1.+D]%*(W-W0)X(W-WO)) 
G72le.-W*W/(Wl1*Wl) 
Pl1=((1.+D1%*(H-W0)%*(W-W0))%*2.*D1+4o%*D1*(W-WO)*DI*(W-W0))/(G2*G2) 
P2=(-(Bo*D1*%* (W-W0)*(W-W0))*D1)/(G2%G2) 

P32 (-(CleD1*(W-WO)*(W-WO))*2,*D1-4.*D1*(W-WO)*D1*(W-WO))/(G2*G2) 
P4= (G1*86*D1*(W-WOD?)*D1*(W-WO))/(G2*G2*G2) 
P9zFI2*(-2./(W1*W1*GT7*GT7)*8. *W*XW/(W1*Wl*W1XW1XG7*GT7*G7)) 
HUZ-XLAM*VX*SINUTXMU*VXCOSU-TG3*G5*(.,25*SINUMQ) 
HUUZ-XLAM*VZCOSU-XMUXVEXSINUTG34G5X(4,25X* COSUMQ) 
HwzG3*(i1s.-.a25*(COSUMQ) )*G4+FI2X*(-20* (W/W1)*(1./W1))/(57%G7) 
HWU=63* (2.2581 SINUMOQ) )*G4 

HwW=  G3%*((l.—-e25*COSUMQO)*(P1+P2+P3+P4))+P9 

DEN-HWU*XHWU -HUU*HWW 
GUM1z-HW*HWU-HUXHWW 

GUM22-HU*HWU-HW*XHUU 

CONTINUE 

RETURN 

END 
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TI 
BS l 


123 
τι 


SUBROUTINE SEARCHIA1sB1sC1,D51,FI1 F12,XLAM O X MÜ 5 COSU ο w I 3 


+WOsG1l9G29W1 »COSUMQ ) 


DIMENSION VINT(10) sWINT(10)sH( 50950) 
FIV=(A1+81*Vx*V)*%(.,01)*(l.-— .25*COSUMO)*G1/G2 
F2V=1./(l.+(W*W/(W1*W1))) 

HAM =XLAM*V*COSU+XMUXVASINU+FIVFFI1+F2VXFI2 
VINT(0)=0.0 

WINT(0)20,0 

DO 131 M=155 

VINT(M)=VINT(M-1)4+3-0 

DO 122 -N=1,5 

WINT(N)=WINT(N-1)+200,0 
G1=C1+D1*(WINT(N)-W0)*(WINT(N)-W0) 
G2=1++D1*(WINT(N)-W0)%(WINT(N)-WO) 
F1IV=(41+B81*VINT(M)*VINT(M))*(.01)*(10-.29*COSUMO )*G1/G2 
rF2Vle/((le*WINTON) *WINTCNO)/(W1*W1)) 
HOMSN)ZXLAM*XVINT((M)*XCOSU*XMU*VINT(M)*SINUTFIV*FIIlITF2VX*FI2 
CONTINUE 

CONTINUE 

AMSRCH=H( 191) 

VSRCH=VINT(1) 

WSRCHZWINT(1) 

DO 121 M=1>5 

DO 123 N=195 

IF (HAMSRCH=H(MsN)) 12331235111 
HAMSRCH=H(MsN) 

VSRCH=VINT(M) 

WSRCH=WINT(N) 

CONTINUE 

CONTINVE l 

IF (HAMSRCH-HAM) 15252 
V=VSRCH 
W=WSRCH 

RETURN 

END 
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SUBROUTINE WORSI (UsVsWsXLAM>XMUsFIl»Al»Bl»C1l>»D1,WO»FI2,»W1»0) 
W=500. 
SINU=SINF(U) 
COSUMQ-COSF(U-Q) 
SINUMQzZSINF (U-Q) 
G1=C1+D1* (W-WO)*(W-WO) 
G2=1.+D1* (W-WO)*(W-W0O) 
G3=(+C1)*FIi*(Al1+B1*V*V) 
G4z(-Gl1*(2. *D1*(W-WO))-4G2*(2e*D1*(W- 10) 317162862) 
G5=(C1+D1* (W-WO)*(W-WO0))/11o+D1%*(W-WO0)*(W-WO) ) 
G7=1le+rW/WO 
PlzCCl.0-D1*(W-WO)* (W-WO))*2,*D14-4. *D1*(W-WO)*D1*(W-WO))/(G2*G2) 
P22(-(8e.e*Di*(W-WO)*(W-WO))X*XD1)/(G2*G2) 
P392 (-(Cl-Dl1*(W-WO)*(W-WO))*2,XD1-4,.*D1*(W-WO)*D1*(W-WO))3/(G2*G2) 
P4z (Gl1*8,.*Dl*(wW-WO)*D1*(W-WO))/(6G2*G2*G2) 
P9-zFI2*((C7-G7T* C,1/WO) -G6*(12/WO))*24. * (1, /WO) )/(€G7*G7*G7)) 
HV =XLAM*COSUtXMURS INU+ (2.01) *F11*(2*B1*V)*G5*(l.-o25*COSUMO) 
HUTZ-XLAM*V*SINUTXMUX*VXCOSU-TG3*G5*(.25*SINUMQ) 
HVV=(.01)*2.*B1*F11*G5*(1l.-.25*COSUMO) 
HVUz-XLAM*XSINUTXMU*COSU-GS5*.01*FI1*2,.*B1*V*(4.25*SINUMQ) 
EE eS o ror cc MEME s: 
-G3*(1.s—-ea25*(COSUMQ) ) *GA-(FI2*((C67*(.,1/WO) ) -G6* (1. /WO) ) )/(G7*G7) 
E G3%*((lo.-.25* COSUMO)*(P1+P2+P3+P4))+P9 
DEN-HVUXHVU-HUU*HVV 
GUM1=-HV*HVU=+HU*XHVV 
GUMZ=HVU*(=HU)+HVX*HUU 
O E 
DU=GUM1/DEN 
DV=GUM2/DEN 
U=U+DU 
V=V+DV 
COSU=COSF (U) 
SINU=SINF(U) 
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COSUMQZCOSF (U-Q) 
SINUMQ-SINF (U-Q) 

Gl1=C1+D1*(W-W0)*(W-WO) 

G221e-Dl1*(W-WO)*(W-WO) 

G3=(+.01)*F11x*(A1+81%*V*%V) 

G4=(-G1*(2.%*D1%* (W-WO0))+G2*(2.*D1*(W-W0)))/(G2*G2) 
G5=(C1+D1*(W-WO)X*(W-W0))/(1.+D1* (W-WO)*(W-W0) ) 
FIV=(A1+B1*VxV)*(.01)%*(l.-.25*COSUMQ)*G1/G2 

F2Vzll. +. 1*W/WO)/(1l. + W/WO) 

HV=XLAM%COSU+XMU*S INU+(.01)%*F1]1%*(2.*B1*V)*G5*(l.-.25*COSUMOQ) 
HUZ-XLAM*XVXSINU-TXMUX*VX*XCOSU-G3*XG5*(.,25*SINUMQ) 
HVV=(.01)*2.*B1*F11*G5%+(1l.-.25*C0SUMOQ) 
HVUZ-XLAM*SINU-XMU*COSU-«-G5*,.01*FI1*2,*B1*V*(,.25*SINUMQ) 
HUUZ-XLAM*VXCOSU-XMU*VXSINU-G3*G5*(4,25*COSUMQ) 
DEN=HVU*HVU-HUU *HV V 

GUM1=-HV*HVU+HU*HVV 

GUM2=HVU* (-HU ) -HVX*HUU 

CONTINUE 

RETURN 

END 
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Appendix II 


In this appendix additional information about the 
paths described in section 7 and the path discussed in 
section 8 will be given. For the first three paths 
the position (x,y) and the controls u,v,w for tho sub- 
WT will be given at twelve-hour intervals. This 
information will be en for Path IV and in addition 
the submarine's location and controls will be included 


for the time steps immediately preceeding and following 


tho corner, 





Path I 

Time x X u X y 
0.0 0.0 οὐ er). 10.7 2076 
19:0 108. LO ου-- ο 0725 
2.0 214,8 = 160% =0.5 10.7 207.8 
36.0 733030 = 197.2 ΠΤ 10.7 207.8 
18.0 449.2 - 244.8 =~ 0.3 20,7 207.7 
60.0 571.5 - 283.1 - 0.3 | 10.7 207.5 
72.0 696,1 ~ 31285 -- 10.7 207.3 
34.0 82243 =- 3332.0 `= 0.1 10.6 207.1 
96.0 949.3 - 315.1  - O.1 ος. 205,8 
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FJ 





108.0 
420.0 
132.0 
144,0 
156.0 
168,0 
180.0 
. 192,0 
204.0 
216.0 
228.0 
240.0 


Time 


12.0 
24.0 
36.0 
49.0 
60.0 


72.0 


τ 


1076.5 
1203.5 
192948 
1555.0 
3598.9 
1701.3 
1822.0 
1941.0 
2058.2 
21738 
2287.2 


2400.0 


1x 


0.0 
102.4 
240.1 
322.4 
438.6 
558.0 
680,0 


276.7 
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Y, 2 
349.2 0.0 
345.7 0.1 
335.1 0.1 
318.0 0.2 
294.8 0.2 
266.0 0.3 
239204 0.3 
193.6 0.3 
150.8 0.4 
104.2 0.4 

- 54.1 0.4 

0.0 0.5 

Path II 
ὦ 2 
040 = 046 
- 66.1 =- 0.5 
124.3  - 0.5 
174.3 = 0.4 
216.3 = 08 
250.4 = 0.2 
- 0.2 


IS 


10.6 
10.6 


10.5 


10.5 


10.5 
10 v5 
10.4 


10.4 . 


10.4 
10.4 
10.3 
119 


IS 


10.1 
10.2 
10.2 
10.3 
10.3 
10.4 
10.4 


ix 


206.5 
206.2 
205.9 
205.5 
205.2 
204.9 
204.6 
204,3 
204.0 
203.7 
203.5 


zB 


1000,0 
1000.0 
1000.0 
1000.0 
1000.0 
1000.0 
1000.0 


em ` 
12388 


96.0 
105,0 
120.0 


13440 


15.0 


156.0 
168.0 
180.0 
192,0 
204.0 
21630 
228.0 


240,0 


0.0 
12,0 
24.0 
36.0 
48.0 


πο 
M 
κ μ.ο 


804.0 
929 4 


1055.7 


1182.4 
1309.2 
1555.7 
1561.7 
1686.9 
1811.2 
1934.5 
2056.6 
2177.6 
2290.2 


2400.0 


0.0 
106.1 
2300 
334.9 
55.8 


X SE 

- 295.6 - 0,1 
- 307-4 = 0.1 
- 312.3 0.0 
- 3109 0.0 
- 303.7 0.1 
- 290.9 0.1 
dou 0.2 
- 250,8 0.2 
> Bald 0.2 
- 195.7 0.3 
- 159.8 0.3 
- 122.7 0.3 
- 64.0 0.5 
0.0 0.5 

Path III 

x u 
0.0  -'0,7 

- 75:4 - 0.6 
~ 144.9 = 0,5 
- 199.1 =- 0. 
- 2456.9 - 0.3 


6/ 


I< 


10.5 
10.5 
10.5 


10.6 


10.6 


10.6 
10.6 
10.6 
10.6 
10.6 
10.6 
10.5 
10.6 


20. 7 


l< 


10.8 
10.9 
10.9 
10.8 


10.8 


a 


1000.0 
1000.0 
1000.0 
1000.0 
1000.0 
1000.0 
1000.0 
1000.0 
1000.0 
1000.0 
1000.0 
1000.0 

705.2 

650.1 


l= 


527.4 
526.4 
526.6 
528.1 
22097 


Timo 


60.0 
72 aD 
84.0 
96.0 
108.0 
120.0 
132.0 
144.0 
159.0 
168.0 
|. 180.0 
192.0 
204.0 
216.0 
228.0 


240.0 


0.0 
- 12,0 


x 


579.6 

705.5 

832.8 

960.6 
1088.3 
1215.4 
1341.5 
1466.3 
1589.5 
1710.9 
1830.5 
1948.1 
2063.8 
2177.6 
2289.6 
2400.1 


[54 


0.0 


70.2 


~ 285.4 
τίς 
- 334.8 
- 356.5 
- 250.1 
- 366,1 
- 335.0 


217.5 
222.9 
264.8 


ο... 
159, 
103.0 


-~ 53.2 


0.0 


ls 


0.3 
0.2 
0.1 
0.1 
0.0 
0.1 
0.1 
0.2 
0.2 
0.3 
0.3 
ο” 
0.4 
0.4 
0.4 


' 0.5 


Fath IV 


Y 


0.0 


- 111.9 
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u 


`- T. 
- 1.0 


I< 


"070 


10.8 
10.7 
10.7 
10.6 
10.6 
10.5 
10.5 
10.4 
10.4 
10.3 
10.3 
100 
10.2 
10.2 


10.4 


I< 


1150 
11,0 


Iz 


534.6 
539.6 
545.9 
553.3 
561.9 
571.7 
582.7 
594.9 
608.3 
622.9 
638.7 
655.8 
67.2 ' 
695.9 
71540 


I< 


200.8 
200,6 





149,6 
237.7 
221.1 


437.9 


5589, 
69,6 
785.7 
910.5 
1038. 3 
1168.3 
1302.0 
1453.5 
1482.9 
1511.5 
1539.4 
1595.1 
1591.5 
1512.3 
1759.6 
1889.7 
2017.5 
2145.1 
220726 
2400.0 


218.0 
217.5 
109.4 
493.1 
568.0 
633.6 
689.6 
736.0 
Veer 
799.8 
816% 
gne 
804.7 
796.0 
784.8 
271.8 
737.2 
742.0 
632.3 
507.3 


581.14 


254.1 


12730) 


0.1 
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0.9 
0.8 
0.7 
0.6 
0.6 
0.5 
0.4 
0.3 
0.2 
0.2 
9, 
0.2 
0.2 
0.2 
0.4 
0.5 
0.5 
0.6 
0.7 
0.8 
0.8 
0.8 
0,8 
0,8 


id 


200.8 


^" 200.8 


200.8 
200,8 
200.9 
200,8 
200,8 
200.8 
200,8 
200,8 
201.1 
202.9 
204.5 
207.0 
211.4 
1000.0 
1000.0 
1000.0 
1000.0 
1000.0 
1000.0 
1000.0 
1000.0 


1000.0 
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